clc,clear;
close all;
%导入相关数据
data=xlsread("random_cities_coordinates.xlsx");
% 城市坐标
cities=data;
% 城市数量
nCities = length(cities);
%计算距离矩阵
dist = zeros(nCities, nCities);
for i = 1:nCities
    for j = 1:nCities
        dist(i, j) = sqrt(sum((cities(i, :) - cities(j, :)).^2));
    end
end
%创建QUBO矩阵
gamma=1000;             %惩罚因子
starcity=5;               %设置初始城市
Q=QUBOTSP(dist,nCities,gamma,starcity);
%退火
[bestx,bestf,bd]=qubosa(dist,nCities, Q,starcity);
disp("Best Solvution");
disp(bestx);
disp("Best Function");
disp(bestf);
tusi(bd,nCities,cities,bestx);


%% QUBO矩阵
function Q=QUBOTSP(distMatrix,numCities,gamma,startCity)
    % 构建 QUBO 矩阵，考虑固定起始和结束城市的约束
    Q = zeros(numCities * numCities);
    
    % 添加路径长度的权重
    for i = 1:numCities
        for j = 1:numCities
            if i ~= j
                % 将距离矩阵的值添加到目标函数中
                Q(i + (j-1) * numCities, j + (i-1) * numCities) = distMatrix(i, j);
            end
        end
    end
    
    % 确保每个城市在路径中只出现一次
    for i = 1:numCities
        for j = 1:numCities
            if i ~= j
                Q(i + (j-1) * numCities, i + (j-1) * numCities) = Q(i + (j-1) * numCities, i + (j-1) * numCities) + gamma;
            end
        end
    end
    
    % 确保路径中没有重复的城市
    for i = 1:numCities
        for j = 1:numCities
            if i ~= j
                Q(i + (j-1) * numCities, j + (i-1) * numCities) = Q(i + (j-1) * numCities, j + (i-1) * numCities) + gamma;
            end
        end
    end
    
    % 固定起始城市和结束城市，并排除在其他位置出现
    for j = 2:numCities
        % 使用一个大值而不是无限大，避免计算出现 NaN
        Q(startCity + (j-1) * numCities, startCity + (j-1) * numCities) = 10^8;
    end
end



%% QUBO退火
function [bestRoute, bestft,bd] = qubosa(dist,nCities, Q,starcity)
    %参数设置
    maxn=1000;       %迭代次数
    T=1000;          %初始温度
    rate=0.99;      %温度变化率
    % 初始化量子态（二进制向量表示路径）
    currentSolution = generates(nCities, starcity,dist);
    currentSolution = enforce_constraints(currentSolution, nCities, starcity);
    % 最优解初始化
    bestx = currentSolution;
    bestf = calcost(currentSolution, Q);
    for iter = 1:maxn
        % 更新量子态（基于退火过程）
        for i = 1:length(currentSolution)
            % 忽略固定的起始和结束城市
            if is_fixed_position(i, nCities, starcity)
                continue;
            end

            deltaEnergy = 2 * currentSolution(i) * sum(Q(i, :) .* currentSolution') - Q(i, i);
            if deltaEnergy < 0 || rand() < exp(-deltaEnergy / T)
                currentSolution(i) = 1 - currentSolution(i);
            end
        end
        % 确保满足约束条件
        currentSolution = enforce_constraints(currentSolution, nCities, starcity);
        
        % 计算新解的代价
        currentCost = calcost(currentSolution, Q);
        
        % 如果是更优的解，则更新最优解
        if currentCost < bestf
            bestx = currentSolution;
            bestf = currentCost;
        end
        % 退火温度逐渐下降
        T = T * rate;

        bestR(iter,:) =biroute(bestx, nCities,starcity);
        bd(iter)=distances(dist,bestR(iter,:));
%         bd(iter)=bestf;
    end
    % 将最优二进制向量解转换为城市路径
    [bestft,index]=min(bd);
    bestRoute=bestR(index,:);
end
function cost = calcost(solution, Q)
    % 计算 QUBO 解的成本
    cost = solution' * Q * solution;
end
function solution = enforce_constraints(solution, numCities, startCity)
    % 确保起始和结束城市固定
    solution(:) = 0; % Reset solution
    solution(startCity) = 1;
    solution(end - (numCities - startCity)) = 1;
    
    % 确保每个城市只出现一次
    for i = 1:numCities
        if sum(solution((i-1) * numCities + 1:i * numCities)) ~= 1
            idx = randi([1, numCities], 1);
            solution((i-1) * numCities + idx) = 1;
        end
    end
end
function fixed = is_fixed_position(index, numCities, startCity)
    % 检查当前索引是否对应于固定的起始或结束城市
    fixed = (index == startCity) || (index == numCities * numCities - (numCities - startCity));
end
function route = biroute(binarySolution, numCities,startCity)
    % 初始化路径，起始城市固定为 startCity
    route = zeros(1, numCities + 1); % 包括返回起始城市的闭环
    route(1) = startCity;
    
    % 用于跟踪哪些城市已被访问
    visited = false(1, numCities);
    visited(startCity) = true; % 标记起始城市已访问
    
    % 逐个确定路径中的城市
    for i = 2:numCities
        found = false;
        for j = 1:numCities
            if binarySolution(j + (i-1) * numCities) == 1 && ~visited(j)
                route(i) = j;
                visited(j) = true;
                found = true;
                break;
            end
        end
        if ~found
            % 如果未找到有效城市，则强制选择一个未访问的城市
            unvisitedCities = find(~visited);
            route(i) = unvisitedCities(1);
            visited(unvisitedCities(1)) = true;
        end
    end
    
    % 最后返回起始城市
    route(numCities + 1) = startCity;

end
% 计算路径的总距离
function distance = distances(dist,path)
    n = length(path);
    distance = 0;
    for i = 1:n-1
        distance = distance + dist(path(i), path(i+1));
    end
%    distance = distance + dist(path(n), path(1)); % 返回起点
end

function initialSolution = generates(numCities, startCity,distMatrix)
    initialSolution = zeros(numCities * numCities, 1);
    route = zeros(1, numCities);
    route(1) = startCity;
    visited = false(1, numCities);
    visited(startCity) = true;
    
    for i = 2:numCities
        currentCity = route(i-1);
        % 找到距离当前城市最近的未访问城市
        distances = distMatrix(currentCity, :);
        distances(visited) = Inf;
        [~, nextCity] = min(distances);
        route(i) = nextCity;
        visited(nextCity) = true;
    end
    
    % 转换为二进制向量表示的初始解
    for i = 1:numCities
        initialSolution(route(i) + (i-1) * numCities) = 1;
    end
end


%% 作图

function ktu=tusi(l,nCities,cities,bP)
l=[l(2:end)];
n=length(l);
count=1:n;
x=zeros(nCities+1,1);
y=zeros(nCities+1,1);
for i=1:nCities
    x(i)=[cities(bP(i),1)];
    y(i)=[cities(bP(i),2)];
end
x(i+1)=cities(bP(1),1);
y(i+1)=cities(bP(1),2);

figure(1)
plot(count,l);
xlabel('迭代次数');
ylabel('最优值');
title('模拟迭代图');
hold on
figure(2)
plot(x,y,'-ro', 'LineWidth', 2, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
xlabel('X 坐标');
ylabel('Y 坐标');
title('路径图');
grid on; % 打开网格
end
